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Abstract 


v  It  has  been  theoretically  shown  that  exponential  rate  of  convergence  in 
the  energy  norm  can  be  obtained  by  a  proper  combination  of  the  h-  and  p- 
extension  of  the  finite  element  method  for  essentially  all  engineering 
problems  in  plane  elasticity,  including  those  with  singularities  in  the  exact 
solution.  In  this  paper  we  present  an  expert  system  frame,  which  guides  the 
user  to  an  optimal  selection  of  the  mesh  and  degree  distribution.  After  a 
crude  preliminary  computation  the  expert  system  predicts  the  performance  of 
various  mesh  and  degree  combinations  and  provides  the  analyst  with  rational 
support  for  his  decisions  about  the  mesh  design.  Numerical  examples  show  that 
the  suggested  approach  gives  superior  results  with  the  minimum  effort  in  human 
data  preparation  and  computer  time. 


Keywords:  Finite  element  methods,  hp-version  of  the  finite  element  method; 
hp-extension,  expert  systems,  elasticity , "'adaptive  methods',  error  estimations. 
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INTRODUCTION 


Today,  powerful  and  sophisticated  methods  in  computational  mechanics  are 
widely  available.  Yet,  these  tools  are  very  often  not  used  as  efficiently  and 
intelligently  as  they  could  be.  The  reason  is  that  the  more  complex  and 
general  this  software  is,  the  more  expertise,  i.e.,  experience  and  knowledge 
of  the  user,  is  necessary  to  use  all  its  possibilities  efficiently.  In 
general,  it  is  impossible  for  a  single  user  to  have  all  this  expertise. 

Traditionally  an  expert's  knowledge  has  been  made  available  through 
direct  consultation  or  by  his  training  of  a  group  of  apprentices.  Both 
approaches  are  slow  and  expensive  and  have  obvious  limitations. 

Today's  technology  of  building  knowledge-based  systems  offers  an  oppor¬ 
tunity  to  change  these  approaches  by  capturing  the  desired  expertise  and 
storing  it  for  use  in  a  computer  system  so  that  newly  acquired  knowledge  can 
easily  be  incorporated  and  efficiently  disseminated. 

The  frame  of  an  expert  system  belongs  to  the  field  of  artificial  intelli¬ 
gence  where  vast  literature  is  available,  (see,  e.g.,  [1],  [2].)  This  frame 
has  been  applied  more  or  less  successfully  in  various  fields.  For  example,  we 
refer  to  [3]  for  a  survey  of  literature  related  to  engineering  problems,  and 
to  [4]  for  the  use  of  an  expert  system  helping  to  select  particular  numerical 
methods.  Various  conferences  and  symposia  have  recently  addressed  this  type 
of  question  in  different  contexts. 

One  of  the  main  parts  of  the  successful  design  of  an  expert  system  is  to 
formulate  (and  justify)  the  rules,  heuristics,  experience  (in  general  exper¬ 
tise)  which  can  be  incorporated  into  the  expert  system  frame.  The  formulation 
of  these  rules  in  the  context  of  an  optimal  mesh  design  for  finite  element 
analysis  will  be  the  topic  of  this  paper.  We  will  also  utilize  the  ideas  of 


adaptive  approaches,  and  we  will  consider  especially  the  problem  of  the 
optimal  use  of  the  possibilities  of  h-,  p-  and  hp-extension  when  the  goal  is 
to  get  the  solution  in  the  range  of  a  prescribed  tolerance  of  the  error  in  the 
energy  norm.  In  a  very  basic  sense,  adaptive  procedures  and  reliability 
assessments  (e.g.,  a-posteriori  estimates)  have  clear  connection  to  the  expert 
system  area  although  they  have  only  algorithmic  character.  (For  a  survey  of 
basic  achievements  in  this  area  and  literature  we  refer  to  [5]). 

Any  finite  element  analysis  can  be  essentially  separated  into  two  parts, 
the  decision  phase  and  the  execution  phase  which  carries  out  the  decisions, 
although  these  phases  are  often  interwoven. 

The  decision  phase  includes  the  specif ication  of  the  goals  of  the  analy¬ 
sis,  (e.g.,  computation  of  stresses  in  some  regions,  computation  of  the  stress 
intensity  factors,  frequencies  of  vibration,  etc.),  specification  of  the 
available  information  and  its  reliability  (as  the  domain  loads,  the  material 
properties,  etc.),  the  desired  accuracy  (as  used  norms  of  the  error,  accept¬ 
able  tolerances)  and  the  available  resources  of  computer  time  and  storage. 
Based  on  these  specifications,  the  user  selects  the  mathematical  model,  con¬ 
structs  the  mesh  using  available  mesh  generators,  selects  the  methodology  of 
computation  of  the  desired  data  (e.g.  elements  and  degrees).  He  selects  the 
method  of  interpreting  the  results  (for  example  graphical  postprocessing)  so 
that  human  and  computer  cost  are  expected  to  be  as  small  as  possible.  How¬ 
ever,  most  finite  element  programs  and  mesh  generators  provide  little  rational 
support  for  proper  decisions  and  obviously  expertise  is  needed  for  an  optimal 
use  of  these  programs. 

There  are  three  basic  versions  of  the  finite  element  method,  the  h- ,  p- 
and  hp-version.  In  the  classical  h-version  the  polynomial  degree  p  is  fixed 
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at  low  levels,  typically  p=1,2,  and  the  accuracy  is  achieved  by  refining  the 
mesh.  The  p-version  fixes  the  mesh  while  achieving  accuracy  by  increasing  the 
degrees  of  the  elements.  The  hp-version  is  a  combination  of  both.  The  theory 
of  the  p-,  and  hp-version  is  given,  for  example,  in  [6],  [7]  [8]  and  [93,  and 
for  practical  applications  we  refer  to  [10]  and  [11].  The  p-  and  hp-versions 
are  recent  developments  and,  besides  some  research  codes,  the  authors  know  of 
only  two  commercial  programs  based  on  the  p-version.  These  are  the  computer 
program  PROBE  (Noetic  Technology,  St.  Louis  ,  USA)  [12]  for  solving  two- 
dimensional  problems  (with  p=1,...,8)  and  FIESTA  (ISMES,  Bergamo,  Italy)  [133 
for  three-dimensional  problems  (p=1  , . . .  ,1)). 

One  of  the  main  features  of  the  hp-version  is  that,  in  practically  all 
engineering  cases  (including  domains  with  corners,  cracks,  etc.)  exponential 
rate  of  convergence  can  be  achieved  (see  [8]). 

Although  theoretically  the  method  is  reasonably  understood  today,  there 
is  still  the  practical  problem  of  how  to  design  the  mesh  and  the  degree  of  the 
elements  so  that  the  given  accuracy  is  obtained  in  the  optimal  way.  The  mesh 
design  and  the  assignment  of  the  degree  of  the  elements  can  be  based  on  an 
adaptive  approach.  We  refer  to  [7]  for  the  analysis  in  a  one-dimensional 
setting. 

The  approach  discussed  in  this  paper  is  essentially  different  from  the 
adaptive  methods,  although  there  are,  necessarily,  similarities  to  the  aims  of 
the  adaptive  codes.  Our  approach  clearly  separates  the  phases  of  decision  and 
execution  and  adds  the  part  of  the  advise  and  communication  with  the  user 
(utilizing  ideas  of  knowledge-based  systems).  It  uses  well-defined  mathemati¬ 
cal  rules,  heuristics  and  experience  to  guide  the  analyst  to  the  design  of 
nearly  optimal  meshes  and  degree  of  the  elements  leading  to  the  desired 
accuracy.  The  main  idea  of  the  approach  presented  here  is  that  after  the 
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design  of  the  most  crude  mesh  modelling  the  geometry  (which  always  has  to  be 
made  by  user's  interference )  one  crude  (cheap)  computation  is  made  to  obtain 
the  main  characteristics  of  the  solution.  Using  this  information,  the  user 
obtains  nearly  optimal  meshes  and  degree  of  elements  (for  the  desired 
accuracy)  together  with  information  about  the  relationship  of  computer  cost 
and  predicted  accuracy. 

Based  on  the  user's  response,  the  nearly  optimal  mesh  and  degree  of 
elements  is  designed  in  one  step  (and  not,  as  in  the  adaptive  codes  mentioned 
above,  in  a  subsequent  mesh  refinement).  The  solution  is  then  obtained  in  the 
execution  phase,  which  is  typically  about  80  to  90J  of  the  total  cost. 

The  knowledge-based  system  ideas  we  are  using  in  our  approach 

a)  help  to  identify  the  class  to  which  the  concrete  problem  belongs, 

b)  characterize  the  critical  areas  of  the  solution  (heuristics), 

c)  give  the  user  concrete  ways  for  extracting  essential  information 
about  the  problem  (algorithms), 

d)  give  the  user  qualitative  information  about  predicted  features  of 
various  optimal  combinations  of  meshes  and  p-distribution  and  the 
relation  between  the  cost  and  accuracy  (heuristics  and  algorithm), 

e)  use  the  analyst's  decision  to  construct  (automatically)  the  mesh,  p- 
distribution,  etc.,  and  present  the  solution  to  the  user  (algorithm) 
with  a-posteriori  assessment  of  the  reliability  (algorithm,  heuris¬ 
tics)  of  the  solution  which  gives  him  the  choice  of  accepting  or 
rejecting  the  result. 

In  this  paper  we  will  restrict  ourselves  to  the  class  of  elasticity 
problems  on  polygonal  domains  with  homogeneous  and  isotropic  material  and 
smooth  boundary  conditions,  although  the  basic  concepts  above  can  be  applied 
to  more  general  situations. 


h 


2.  THE  ELASTICITY  PROBLEM  ON  A  POLYGONAL  DOMAIN  AND  ITS  BASIC  PROPERTIES 


Let  us  consider  the  plune  (strain)  elasticity  problem  on  a  polygonal 


domain  ft  shown  in  Fig.  2.1. 


4  *3, 


^2  v  ^2 


I  (A, 


Figure  2.1.  Polygonal  domain. 


By  3ft  we  denote  the  boundary  of  ft  and  have 

3ft  =  J  r.;  Y  *  { 1 , . . . , M }  where  r.  are  the  open  edges  connecting  the 

i £  Y  1 

vertices  A^f  A^+1  (AM+.,  =  A-|  ) .  By  u>  we  denote  the  internal  angles  of 

in  A;  .  We  will  not  exclude  the  case  oj.  =  2ir  and  w.  =  tt  . 

“  l  l 

Let  3ft  »  rD  'J  rN  ; 

fD  “  i  D  c  Y  Fi  ;  rN  =  "  ~D  ' 


TTTT1 


Fd  will  be  called  displacement  and  F^  traction  boundary.  We  allow 


rD  ■  3o 


or  rc  - 


Considering  the  plane  (strain)  elasticity  problem  on  P.  ,  we  denote  by 
u— ( u ^ ,  U2)T  the  displacement  vector,  E,  0  £  v  <  1 /2  the  Young's  elasticity 
module  and  Poisson's  ratio,  respectively.  The  strain  energy  is  given  by 


9u,  2 


W(u)  =  2( 1-v) (1 +v)  l  ^(1  v)  ^9x^  +  2v  3x1 


9u,  9u, 


9x, 


(2.1  ) 


9u  2  9u  9u  2 

(l-v)  (^-)  *  (-P)  *  5^)  ]  IX 


on  Fd}  be  the  energy  space  and  let  |u|‘ 


Let  E  —  {U  |  W(U  )  ^  00  ,  U  —  0  Ull  i  |^J  anu  -.v--.  |  -  | 

W(u).  It  is  well  known  that,  the  energy  space  E  is  equivalent  to  the  Sobolev 

space  ( P)  (modulo  rigid  body  motion  displacement  functions  when  rD  =  0  ). 

The  elasticity  problem  (without  volume  forces)  for  prescribed  tractions 
T  =  (t^ ,  t2)T  on  rN  now  has  unique  (weak)  solution  u0  which  minimizes 
over  E  the  quadratic  functional 


(2.2) 


n(u)  =  W(u)  - 


7  f  t . u . ds 
*•  J  ii 

i=i  r. 


N 


If  fp  =  0  then  we  assume  that  the  tractions  satisfy  the  usual  equilibrium 


condition.  Furthermore,  we  assume  that  t^  i=1,2  are  analytic  on  every 


(closed)  edge  F^  c  .  The  solution  Uq  is  then  analytic  on  7?  -  ^  A ^ 

and  for  some  0<B<1  ,  d>1 ,  c>0,  I  a | ^  2  . 


(2.3) 


[  j  (Dau0>i)2[4>]2(!ai‘2+6)dP]1/2  £  Cd  f a  i  |  c  |  ! 


-  ■- 


-  V-  V  -^-e  ^ 


holds  for  all 


a  =  (a1 ,a2)  ,  where  a  £  0  integer,  =  |a| 


$  =  n  jx-A.J.  | x — A .  |{  is  the  Euclidean  distance  of  x  to  A^,  and 

i-1  1  1 


Dau  = 


<V  a2. 


a-i  ct2 


.  Further,  in  the  neighborhood 


3  Xi  3  x^ 


S  =  Q  |x  |  |x-A.  | <r Q }  of  the  vertices  Ai  the  exact  solution  uQ  can  be 


written  in  the  form 


(2.4) 


l  K..4>..(r.)g..(0.)  +  w  . 

A  ij  ij  i  ij  i  Qj.i 


where  (r ^ ,  0.)  are  polar  coordinates  with  respect  to  Ai ,  <J>. .(r.)f  g  (0.  )  are 
a-priori  known  functions  independent  of  uQ  (depending  on  uk)  and 


(scalars)  are  stress  intensity  factors  (dependent  on  u0  )  and  w 


Q.,  i 


is  an 


analytic  function  on  S.  -  A.  which  is  smoother  than  the  first  term  on  the 

ii 


right-hand  side  of  (2.4).  is  a  positive  integer  which  will  be  chosen  for 


our  practical  purposes  as  1  or  2.  For  more  details,  see  [I1*],  [15]  and  [16]. 


function  g. .(0^)  is  analytic  up  to  the  boundary  of  Si  and 


^j(r^)  =  Re  (ri  1")2-gp  r  ^ )  ;  p  2  0  integer.  ok  .  are  in  general  complex 


numbers  with  positive  real  parts  and  Re  a.  .  ,  5  Re  .  .  Hence  the  first 

F  i,J  +  1  ij 


terms  in  (2.4)  are  the  most  singular  ones.  For  the  way  to  compute  stress 


intensity  factors  we  refer  to  [17],  and  [ 1 8 ] . 


Table  2.1  shows  the  values  of  a.,  and  a._  for  some  internal  angles 

ii  1 2 

uk  when  the  tractions  are  prescribed  on  both  sides  of  A  In  this  case 


p  =  0 
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3.  THE  HP-VERSION  AND  ITS  BASIC  PROPERTIES 


As  was  stated  in  the  introduction,  the  hp-version  is  a  proper  combination 
of  the  h-  and  p-version  of  the  finite  element  method.  For  simplicity  we  will 
assume  that  the  polynomial  degree  of  all  the  elements  is  the  same.  It  has 
been  shown  [8]  that-  if  the  exact  solution  satisfies  (2.3)  then  for  a  proper 
mesh  A(p)  depending  on  p  the  error  of  the  finite  element  solution  uEE 
measured  in  the  energy  norm  decays  exponentially.  More  precisely 

(3.D  |e|E  =  |uFE  -  uQ|E  <  Ce"a  N  (p;  P  *  * 

where  a>0  and  N  is  the  number  of  degrees  of  freedom. 

The  error  in  particular  elements  has  two  essential  parts.  The  local  one 
depends  only  on  the  solution  in  the  given  element  and  the  global  one  reflects 
the  influence  of  the  error  in  the  entire  domain  (so-called  pollution  error). 
The  global  error  is  (for  properly  designed  meshes)  smaller  than  the  local  one. 
If  an  element  is  not  adjacent  to  any  vertex,  then  the  local  error  of  the  (p- 
version)  finite  element  method  is  of  the  order  e-a  ^  where  is  the 

number  of  degrees  of  freedom  associated  with  the  element  whereas  the  error  is 

_  O 

of  order  when  the  particular  element  has  a  vertex  in  a  corner  of  the 

domain.  Hence  the  decay  of  the  error  is  more  rapid  in  areas  far  from  corners 
of  the  domain  than  in  the  neighborhood  of  the  corners.  Therefore  only  proper 
refinement  in  the  neighborhood  of  corners  is  needed,  while  in  the  elements  not 
adjacent  to  corners  high  accuracy  can  be  achieved  only  by  increasing  p.  For 
practical  engineering  accuracies  we  can  treat  convex  corners  as  no  corners 
(provided  that  the  boundary  condition  is  the  same  on  both  sides  of  the 


vertex) . 


So  the  proper  mesh  A(p)  is  a  geometric  one  with  geometric  refinement  in 
the  neighborhood  of  every  reentrant  corner.  An  example  of  such  a  mesh  with  2 
layers  at  Aq  is  shown  in  Fig.  3.1  • 


Figure  3.1.  Example  of  a  geometric  mesh  with  two  layers. 

It  has  been  shown  that  the  optimal  number  n  of  layers  of  elements 
increases  with  the  polynomial  degree  p,  that  the  optimal  ratio  of  the 
geometric  mesh  is  independent  of  the  strength  of  the  singularity,  p,  and  the 
number  n  of  layers  and  has  a  magnitude  of  .15.  (see  [7]  and  [8].)  Thus  a 
very  strong  grading  towards  the  singularity  is  obtained. 

If  the  mesh  is  fixed  (with  different  numbers  of  layers)  and  the  polynom¬ 
ial  degree  p  increases,  then  the  error  behaves  as  schematically  shown  in 


Number  of  degrees  of  freedom  N 


Figure  3.2.  jeL  in  dependence  of  the  number  of  degrees 
of  freedom  for  different  number  of  layers  n  . 

For  each  fixed  number  of  layers  we  see  a  typical  reverted  S-curved  be¬ 
haviour  of  the  error.  We  see  the  preasymptotic  phase,  (curved  down),  when  rhe 
error  decays  exponentially  and  the  asymptotic  phase,  (straight  line),  when  the 
error  decays  algebraically.  For  theoretical  results  we  refer  to  [7]  and  [81- 
Figure  3.2  shows  that  at  reentrant  corners  for  high  number  of  layers  and 
low  polynomial  degree  p  the  mesh  can  be  overrefined  (with  respect  to  N  ) 
while  an  insufficient  number  of  layers  leads  to  an  underrefined  mesh.  The 


detailed  behaviour  is  complicated  but  it  strongly  depends  on  the  values  of 

K<  4  and  a.  .  in  (2. 4 ) . 
iJ 

For  practical  reasons  we  wish  to  design  a  mesh  (i.e.  the  number  of  layers 
in  every  corner)  depending  on  p  so  that  the  desired  accuracy  can  be  achieved 
in  a  most  economic  way. 

The  definition  of  the  cost  of  a  computation  can  vary.  It  can,  for  ex¬ 
ample,  be  a  function  of  the  degree  of  freedom  N  or  the  number  of  operations 
required  in  the  computation  and  elimination  of  the  stiffness  matrix,  depending 
on  the  sparsity  of  the  matrix. 

Our  experience  shows  that  (because  of  overhead,  10)  it  is  enough  to 
consider  the  cost  as  a  function  of  N  only. 


4.  THE  MAIN  PRINCIPLES  OF  THE  MESH  DESIGN 


The  main  goals  for  the  optimal  mesh  design  are  the  following: 

1)  The  user  supplies  only  the  basic  data  which  have  to  be  kept  to  a 
minimum. 

2)  The  user  should  get  rational  support  for  his  finite  element  analysis, 
mesh  design,  relationship  between  accuracy  and  cost  of  the  computation,  etc., 
and  he  should  be  given  the  opportunity  to  make  the  crucial  decisions  in  the 
entire  process.  The  expert  system  provides,  then,  the  support  for  the  user  to 
execute  his  decisions  automatically. 

The  finite  element  analysis  based  on  the  expert  system's  support  has  the 
following  components: 

a )  Interactive  graphical  mesh  generation  and  data  characterization 

The  finite  element  process  starts  with  the  characterization  of  the  prob¬ 
lem  by  definition  of  the  domain  by  a  basic  mesh  which  is  as  coarse  as  possible 
(or  by  solid  modeller  information)  of  the  boundary  conditions,  the  loads,  ma¬ 
terial  properties,  etc.  Graphical  input  should  be  used  wherever  possible  and 
analysis  of  the  data  for  possible  input  errors  should  be  provided. 

In  this  paper  we  will  concentrate  only  on  the  aspects  of  the  domain  defi¬ 
nition  and  we  will  assume  that  they  belong  to  the  class  of  polygons.  The 
basic  mesh  will  also  serve  as  the  basic  set  of  elements  which  will  be  dealt 
with. 

b )  Basic  decisions  about  critical  and  noncritical  elements  of  the  basic  mesh 
An  element  is  said  to  be  noncritical  if  the  optimal  rate  of  convergence 

(in  the  range  of  engineering  accuracy)  can  be  obtained  by  increasing  the 


degree  of  the  elements  without  any  refinement.  An  element  is  said  to  be 
critical  if  the  optimal  rate  of  convergence  cannot  be  achieved  only  by 
increasing  the  degree,  but  the  element  has  to  be  refined.  The  decision  about 
these  two  classes  has  to  be  based  on  the  expertise  which  uses  theoretical 
results,  heuristics  and  experience. 

For  example,  all  elements  adjacent  to  reentrant  corners  or  vertices  with 
changing  boundary  conditions  as  well  as  elements  with  excessive  aspect  ratio 
or  elements  where  the  load  is  unsmooth  are  critical.  We  will  assume  in  this 
paper  that  the  critical  elements  of  the  basic  mesh  are  only  those  adjacent  to 
vertices  with  concave  internal  angles. 

c)  Extraction  of  critical  data  about  the  solution 

Crude  preliminary  computation  and  various  extraction  procedures  give  the 
basic  qualitative  and  quantitative  information  about  the  solution  in  critical 
areas.  Especially  in  this  paper,  the  approximate  determination  of  the  stress 
intensity  factors  in  (2.4)  will  give  essential  information  about  the  solution. 
We  will  extract  the  stress  intensity  factors  by  the  procedures  described  in 
[17]  and  [18]. 


d )  Prediction  of  the  performance  of  various  optimal  mesh  and  degree  combina¬ 
tions 

This  prediction  is  one  of  the  critical  parts  of  the  system  and  will  be 
elaborated  on  in  detail  in  section  5.  The  prediction  characteri zes  the  per¬ 
formance  of  the  meshes  with  different  numbers  of  refinement  layers  n  in  the 


critical  elements  in  dependence  of  the  polynomial  degree  p.  This  prediction 


utilizes  the  information  obtained  by  the  extraction  technique  under  c)  above. 


Based  on  the  analysis  in  d),  the  basic  data  about  the  predicted  relation 
between  cost  and  accuracy  are  graphically  presented  to  the  user.  Other  char¬ 
acterizations  may  also  be  provided.  For  example,  in  the  case  of  multiple 
loads  the  user  should  get  advice  on  whether  it  is  advantageous  to  solve  the 
problem  with  one  mesh  or  to  use  different  meshes  for  various  classes  of 
loads.  In  this  paper  we  will  restrict  ourselves  to  the  case  of  one  load  only. 


f )  Mesh  generation  for  the  (final)  computation  based  on  the  user's  decision 
Based  on  the  user’s  decision,  the  mesh  which  is  expected  to  give  accept¬ 
ably  accurate  results  in  an  optimal  way  (for  given  optimal  degree  of  the 


elements)  is  constructed  automatically. 


;)  Finite  element  computation  and  assessment  of  the  reliability  of  the 


computed  data 


The  finite  element  code  uses  the  optimal  mesh  and  determined  degree  of 


the  elements.  It  provides  the  solution  and  a-posteriori  assessment  of  the 
reliability  of  the  basic  computed  data  (e.g.,  error  in  energy  norm).  In  this 
paper  we  will  use  the  code  PROBE  for  the  computation. 


h)  Presentation  of  the  basic  results  of  the  computation  for  the  user's 


acceptance 


The  basic  results  obtained  under  g)  are  presented  to  the  user  and  with 


the  possible  assistance  of  the  expert  system  the  user  decides  whether  the 


results  are  acceptable. 


If  the  user  accepts  the  finite  element  results,  various  postprocessing 
procedures  are  appended  to  give  him  the  results  in  which  he  is  interested,  to 
provide  graphical  display,  etc.  The  expert  system  will  probably  help  to 
interpret  the  results. 

j )  Processing  of  rejected  results 

If  the  user  rejects  the  results  because  the  prediction  was  not  accurate 
enough,  then  the  computed  data  are  used  analogously  to  the  crude  mesh  data  - 
with  additional  information  and  the  process  is  repeated  from  c).  Rejection 
should  occur  only  under  very  exceptional  circumstances. 

The  underlying  data  structure  of  the  whole  process  is  a  project  data  base 
where  the  complete  history  of  the  analysis  is  stored.  This  makes  it  easily 
possible  to  add  new  load  cases  or  to  modify  the  needs  of  the  accuracy  of  the 


results. 


5.  PREDICTION  OF  THE  PERFORMANCE  OF  VARIOUS  OPTIMAL  MESHES  AND  DEGREES  FOR 
CRITICAL  ELEMENTS 


5 . 1  Error  analysis  in  critical  elements 

We  have  seen  that  in  the  neighborhood  S^^  of  the  vertex  Ai  the  solu¬ 
tion  has  the  form  (2.4).  Because  the  function  g,  S(Q)  is  smooth,  the  be- 

f  J 

haviour  of  the  error  is  qualitatively  the  same  as  if  g.-  ^(0)  =  1.  Hence  we 

J-  >J 

can  expec..  that  the  error  in  the  sector  S^  is  essentially  the  same  as  in  the 
one-dimensional  setting  in  the  interval  (0,  rg)  with  the  weight  x  express¬ 
ing  polar  coordinates.  For  the  two-dimensional  analysis  we  refer  to  [8]  which 
justifies  our  assumptions. 

Let  us  first  study  the  problem  of  the  best  approximation  on  the  interval 
I-(-1  ,1 )  .  Let  for  E  <  1 


(x  -  E)  + 


|  x  -  C  for  x  >  E 
1  0  for  x  S  E 


and  let  i|>(a,E,x)  =  (x  -  E). 


Denote 


(5.1  ) 


2  i  1/2 


E  ( a ,  E )  =  [inf  /  (<j/  -  id)  dx] 

P  u  -1 


where  the  infimum  is  taken  over  the  set  of  all  polynomials  of  degree  £  p  . 
Analagously  let 


(5.2) 


Ep  (a,E) 


2  .  il  /2 


[  inf  J  (1  -  x2)(i|;  -  u>)  dx] 

a)  -1 


The  values  Ep(a,E)  and  Ep  (a,E)  express  the  minimal  error  achievable  by 


the  polynomial  approximation.  Now  we  have 


j-  *  ^7  ,iK7.v71p  y^y*y9y?ymywymyw 


**■*  ***  4>~  IW  $tm  t>«  UI  «jVj^J 


Theorem  1 


a)  If  g  =  -1  then 


(5.3) 


where 


Ep(a, 


C0(a) 


?ct  +  % 

"  =  co(a)  — 7^ 

(p+1  ) 


-  ('  *  »<?)) 


(T(1+a))  1  sin  ( Tig) 

it  /  2a+1 


b)  If  5  <-1  then 


(5.M 


where 


Epu.u  -  c,(a)(l^)“-^r  (1  *0  (1^) 

H  (p+1)  p 


C  >  o  ,  r  = 


5  +  /  t  -1 


(a)  = 


r(1+a)  |  sln(ira) 


T(a)  is  the  gamma  function  and  O(-)  is  a  function  which  is  of  the  same 

order  as  ^  .  For  the  proof,  see  [7]  theorems  5  and  6. 

We  see  that  for  E  =  -1  ,  i.e.,  when  the  singularity  is  in  the  vertex, 

the  error  decays  algebraically  with  the  order  (p  +  1)2a+1  ,  while  when  E,  <-1  , 

i.e.,  the  singularity  is  outside  the  interval,  the  decay  is  exponential. 

* 

A  similar  theorem  holds  for  : 

Theorem  2 : 

a)  If  E  =  -1  then 


(5.5) 


where 


Ep  (a,  -1) 


C2(a) 


(p*1  )‘ 


-2(i  •  o(h) 
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k* 

l 

w 


„  ,  s  (r(a+2))  sin(ira) 

C  (a)  =  - V/p - -J- 

2  n(a*1  Y/d 


b)  If  5  <-1  then 


(5.6) 


*  2  r  1 — 01  r^+^  r  d  1 1 

E_  (o.t)  =  C  (a)(1-r2)[^-)  -~  (l+o(--)) 

p  i  2r  (p+l)a  p5 


where  z,  >  0  ,  r  is  defined  as  in  theorem  1  and 


_  ,  v  1  T(a+2)  sin(ita) 

C3 a  ■  - 7= — 


Theorem  2  can  be  easily  obtained  by  modifying  the  proof  of  theorem  1 , 
using  the  fact  that 


.♦1  ,  r  0  for  j  *  k 

( 1  -x  )  P!PV  dx  -  i  ,+1 

h  J  k  l  j(J+1)  J  Pj‘ 


dx  for  j  =  k 


where  P,  are  the  Legendre  polynomials. 

J 


* 

Comparing  (5.3)  and  (5.5)  we  see  that  is  smaller  than  Ep  and  that 

# 

E  0  for  a  >  -1  while  E  -*  0  only  for  a  >  -  1 /2  . 

P  P 

* 

On  the  other  hand,  comparing  (5.^)  and  (5.6)  we  see  that  E^  and  E^ 

* 

are  up  to  a  constant  essentially  the  same,  (i.e.,  Ep  *  E^  ). 

Let  us  now  define 

*#  r  r+^  2  n  i  /? 

E  (a,  O  =  inf  [  j  (1+x)(^-oj)  dx  J 

^  03  ~  1 

*  *  ** 

As  Ep  *  Ep  for  E  <  -1  ,  it  is  clear  that  in  this  case  also  Ep  s  Ep 

#  ## 

Moreover,  Ep  *  Ep  also  in  the  case  E  =  -1  ,  since  near  the 

p 

singularity  -  1  of  ij)  the  weights  1  +x  and  1  -  x  are  essentially  the 
same . 


1? 


So  far  we  have  beer,  interested  in  the  case  I  =  (-1f1).  Let  us  now  take 
I  =  (a,b)  and 


w  * 


d  <  a 


Defining 

E  (a.b.d.a)  »  [  inf  f  ( ip  -  w)2  dx  ]  ' /2 

P 

K  ui  a 

and 

E  *(a,b,d,a)  =  [  inf  {  (x-a )  (b-x)  ( oo ) 2dx  ]  1/2 

P 

y  oj  a 

^  ^  ^  p  4  /  p 

E  ( 0 ,  b  ,  0 ,  a)  =  [  inf  J  x(i|>  -  u)  dx] 

p  u  0 


we  easily  get  the  relationship  between  Ep(a,£)  and  Ep(a,b,d,a)  by  a 
scaling  argument. 

Lemma  3 : 


(5.7a) 

Ep(a,b,d,a)  = 

h“  +  1  /2  Ep  (  ^ ,  a ) 

(5.7b) 

* 

Ep  ( a ,  b  ,  d ,  a ) 

-  h“  +  3/2Ep*(t,a) 

(5.7c) 

E  ( 0  ,b ,  0 ,  a) 

P 

n  +  1  *# 

=  h  E  ( 0 ,  a) 

P 

where 

xTH 

It 

O 

a+b  ,  b-a 

=  —  h  =  — 

Now  let 

A  be  the  one-dimensional  mesh  on  I=(0,B). 

V 

o 

X 

II 

o 

< 

x,  <  x~  <  ...  <  x  *  B 
1  c  M 

where 

x  =  BqM_i 

,  q  <  1 

and  denote  I.  =  (x.  ,,  x.) 

J  J-1  J 

=  x  and  define 


Assume  that  f 


(5.8) 


E(6,p,q,M)  =  (inf  j  (f'  -  oi')2  x  dx)1 

(l)  0 


where  the  infimum  is  taken  over  all  functions  id  (  H'(I)  and  w  is  a 
polynomial  of  degree  p  on  1^,  i  =  and  w(x^)  =  f(x^). 

We  now  have 
Theorem  3 : 


(5.9)  E(B,p,q,M)  =  C-B*  [■ 


28(«-n  r2p  (,  26(H-1))(1  6-,  2 


1  -  q 


where  D.,(6)  £  C  £  D2(6)  with  independent  of  M  and  p  . 

Proof : 

We  will  approximate  f'  =  Bx  by  polynomials  of  degree  (p-1 )  on 
every  Ij  .  Then  by  integration  we  construct  w  -  h\i). 


On  I j ,  j  >  1  we  use  lemma  3  and  set 


a  =  x .  , 
J-1 


BqM+1"J  ;  b  =  x^.  =  BqM  j 


h  =  BqM_j  ^ 


BqM"j  l|a 


5  "  -h 


( 1  +q )  +2/ q 


We  get 


E(x  .  ,x  , 0,  B~1  )  =  C(B)  -  —  —  B(6'V2)q(M_j)(6'1/2)(^)B_1/2 


J-1  ’  J 


(Er)^  p8 


=  C ( 8 )  — r-  ( 1  - q ) *^q  2  '-r  q (M~ j  }  (  8  ‘4) 

O  '2  _  tO 


On  I-,  we  have 
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We  will  assume  that  the  solution  to  be  approximated  is 


K1  r  1  4>1  ( ©)  +  K 2r  24>2(0) 

and  that  the  critical  element  is  divided  into  M-1  layers  by  a  geometrical 
mesh.  Furthermore,  let 

2  ^  p  p  a . 

(5.10)  p  =  C-u>  l  K.^E  (o  ,p,q,M)-B. 

i=1 

to.  is  the  angle  of  the  critical  element  at  Ai  (see  Fig.  (5.1).)  and  B.^  is 
the  radius  of  the  smallest  circle  centered  at  Aj  covering  the  element. 

The  expression  (5.10)  is  an  error  predictor  for  the  error  in  the  energy 
norm.  It  will  be  seen  in  the  numerical  examples  that  p  reliably  gives  all 
characteristics  of  the  error  behaviour  in  dependence  of  M  ,  i.e.,  the  number 
of  layers. 

The  value  q  will  be  chosen  of  the  order  .15  as  stated  above,  which  is 
the  optimal  value  for  one-dimensional  problems.  C(B)  in  (5.10)  is  not 
critical  for  our  purposes  because  it  ranges  only  over  a  small  interval. 

5.2  Optimization  of  the  number  of  layers  in  the  critical  elements 

As  has  been  said,  we  will  assume  that  the  decisive  areas  for  accuracy  are 
the  critical  elements  which  have  to  be  refined.  Hence  we  will  assume  that  the 
only  error  is  in  these  elements  and  that  the  stress  intensity  factors  are 
known . 

For  given  p  and  number  of  layers  n^  =  M^-1  in  the  t(i)  elements 

adjacent  to  the  critical  vertices  A^  ,  i=1 . s  we  get  from  (5.10)  the 

error  predictor  for  the  error  in  the  whole  domain  P 
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(5.11)  p  (p,M  ) 


s  t( i )  2  ?  ? 

>  I  I  ,  I  K  V(o  .p.q,  M  )B  *’ 

i=1  j=i  X*J  £*i  fc>1  1  1>J 


fi.»i 


where  ^  are  the  (known)  stress  intensity  factors  in  the  vertex  and 


otf  .  the  exponents  of  the  singular  functions  in  (2.M). 

*>  i  i 


The  number  of  degrees  of  freedom  over  all  the  refined  mesh  will  be  de¬ 
noted  by 


N(p,M, 


i-1 , . . . , s) 


Assuming  that  the  computational  work  is  a  function  of  N  only,  we  can  formu¬ 
late  the  following  optimization  problem: 

Given  an  upper  bound  Nq  of  the  total  number  of  degrees  of  freedom  , 
find  p  and  so  that  it  minimizes  (5.11)  under  the  constraint 


N(p,M.  ,  i-1 . s)  <  NQ  . 


For  practical  purposes,  p  can  be  restricted  to  a  maximum  degree,  Pmax»  and 


it  is  also  reasonable  to  restrict  the  number  of  layers  to  some  Mmax  .  A 


possible  choice  for  the  maximally  allowed  number  of  layers  is  M_.v  =  2p 

Ilia  A  Ilia  A 


as  it  was  shown  in  [8]  that  the  optimal  combination  of  M  and  p  in  the  case 
of  a  single  crack-tip  singularity  is  given  by  M  =  p  . 

With  these  restrictions  the  optimization  problem  is  finite  and  various 
methods  for  constrained,  finite  optimization  could  be  used.  Yet  it  turns  out 
that  (5.11)  often  has  many  local  extrema  so  that  the  exact  optimum  will  be 
hard  to  find.  Nevertheless  it  is  not  necessary  to  find  the  minimum  exactly 
because  various  simplifications  have  been  made  in  obtaining  (5.11).  In  all 
our  test  examples,  the  following  simple  heuristic  algorithm  gives  nearly 
optimal  results  with  only  minor  computational  effort. 


2h 


Perform  steps  1  to  5  for  p=1 , . . . .Pmax: 

Step  1:  Choose  as  starting  distribution  =  Mmax ,  i=1,...,s 
Step  2:  Compute  p(p,PL  ,  i  =  1,...,s)  and 

N(p ,NL  ,  1  ~  1  ,  . . • , s )  . 

Perform  steps  3  and  4  until  N  £  Nq  . 

Step  3:  Find  j  ?  {1,...,s}  with  the  following  property: 

If  one  layer  is  removed  at  Aj  ,  the  quotient  of  the  change  of  the  total 
number  of  degrees  of  freedom  and  the  change  in  the  error  predictor  (5.11)  is 
maximal . 

Step  4 :  Remove  one  layer  at  Aj  ,  update  pCp.Nh  ,  i*1,...s)  and 

N(p,Mj  ,  i=1 , . . . ,s) 

Step  5:  { ,  i=1,...,s}  is  now  a  locally  optimal  distribution  for  the 

specified  p-degree. 

In  practice  it  is  often  useful  to  compute  the  error  predictor  as  a 
function  of  Nq,  for  optimal  distributions  of  the  degree  p  and  the  number  of 
layers,  i.e.,  to  perform  a  sequence  of  optimization  processes. 

5.3  Prediction  of  the  total  error 

Because  of  the  constant  C  in  (5.10)  the  error  predictor  estimates  the 
total  error  in  the  energy  norm  only  up  to  a  constant,  which  has  still  to  be 

determined.  As  has  been  shown  in  [8],  the  error  behaviour  for  optimal  meshes 

1  /3 

is  of  the  type  je|E  -  e-^N  .  Hence  it  is  advantageous  to  visualize  the 
error  in  the  scale  lg|e(  x  N1  ^  #  The  error  predictor  gives  good  qualita¬ 

tive  characterization  of  the  total  error  (including  the  errors  in  noncritical 
elements),  if  it  will  be  properly  scaled  by  a  shift  of  lg(|e[  ).  This 
adjustment  determines  the  proper  constant  C  in  (5.10)  and  will  be  called 


calibration.  The  main  idea  is  to  calibrate  the  error  predictor  (i.e.,  the 

curve  lg(|e|)  x  N1 ^  so  that  for  small  p  the  error  predictor  is  close  to 
the  actual  error.  This  can  be  done  by  various  approaches.  For  example,  we 
can  compute  the  solutions  for  p=1  and  p=2  or  p=1,2,3  on  the  basic  mesh 
and  determine  the  approximate  error  by  extrapolation,  assuming  that 

(5.12)  |  e  1 E  =  C  N-6 

for  the  p-version  on  a  fixed  mesh  where  6  depends  on  the  smoothness  of  the 
solution.  Experience  shows  that  for  low  p  and  coarse  meshes  6  =  .70  to 
6  =  .75  gives  reasonable  results  in  most  practical  cases. 


6.  NUMERICAL  EXAMPLES 


In  two  numerical  examples  we  will  show  the  steps  of  the  expert  system  as 
described  in  section  4.  The  first  example  will  have  only  one  singular  point, 
showing  the  basic  principles  of  the  error  prediction.  The  second  example  will 
have  9  singular  points  and  show  that  superior  accuracy  can  be  obtained  in 
practical  problems. 

6.1.  Problem  1 ;  The  L-shaped  Domain 

Figure  6.1  shows  the  basic  mesh  of  three  elements  constructed  in  step  a) 
of  section  4.  The  domain  is  loaded  by  tractions  so  that  the  exact  solution 
shows  mode  1  stress  distribution  given,  for  example,  in  [18],  The  solution 
has  singular  behaviour  at  the  reentrant  corner. 

As  each  of  the  three  elements  is  adjacent  to  the  reentrant  corner,  In 
step  b)  each  element  is  expected  to  be  critical. 

Step  c)  extracts  critical  data  about  the  solution  of  the  problem.  In 
order  to  get  sufficiently  accurate  data  for  an  extrapolation  of  the  exact 
energy  (see  section  5.3)  and  an  accurate  relation  of  the  two  stress  intensity 
factors  and  K2  at  the  reentrant  corner,  one  circular  refinement  layer 

was  added  at  the  reentrant  corner  (Fig.  (6.2)).  On  this  mesh  a  finite  element 
computation  was  performed  for  p=1,2,3.  This  crude  computation  yielded  for 
p=2  approximate  stress  intensity  factors  Kn  =  0.91 ,  K2  *  0  (the  exact 
values  for  this  example  are:  K-j  =  1  and  Kj  =0).  The  approximate  energy  of 
the  three  solutions  was  U1  =  3-8288279,  U2  =  4.0975754,  U3  =  4.134477.  The 
stress  intensity  exponents  are  only  geometrical  data  and  can  be  computed 
exactly.  They  are  given  by  o1  =  0.544484  and  a2  =  0.908529  . 
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Number  of  degrees  of  freedom  N 


Figure  6.3.  L-domain.  Predicted  error 

Steps  d)  and  e)  predict  the  performance  of  various  optimal  meshed  and  p- 
degrees  and  present  this  prediction  graphically  to  the  user.  Out  of  the 
computed  energy  of  the  solutions  for  p«1,2,  the  exact  energy  is  extrapolated 

a 

by  (5.12)  with  <5  =  .75  to  U  «  4.16^6802  which  estimates  well  the  exact 
strain  energy  of  Uex  =  4.1 5M 5U42  * U  served  to  calibrate  the  predicted  error  in 
the  energy  norm.  Figure  (6.3)  shows  the  error  prediction  as  it  is  presented 
to  the  user  for  p  =  1 . 8  and  n-2,...,7. 


2? 


In  steps  f)  and  g)  one  particular  mesh-degree  combination  would  be  con¬ 
structed  automatically  after  the  user’s  decision  and  the  finite  element 
computation  would  be  performed.  In  order  to  show  the  performance  of  the  error 
prediction,  we  compare  Fig.  (6.3)  with  the  actually  computed  solutions  on  the 
combinations  of  mesh  and  polynomial  degree.  The  numerical  results  are  taken 
from  [19].  Figure  (6.4)  gives  the  real  error  in  the  energy  norm.  All  mesh- 
degree  combinations,  which  were  predicted  to  be  optimal  for  a  specific  number 
of  degrees  of  freedom,  are  marked  with  circles  in  Fig.  (6.4).  it  shows  that 
each  of  these  meshes  is  on  the  lower  left  envelope,  which  means  that  the 
prediction  of  the  optimality  was  right  in  each  case. 


200  800  1500 

Number  of  degrees  of  freedom  N 


Figure  6.4.  L-domain.  Exact  error. 


The  slope  of  the  curve  joining  the  optimal  mesh-degree  combinations, 
i.e.,  the  envelope  in  Fig.  6.4  is  a  straight  line  in  the  Log-N1  ^-scale ,  which 
shows  that  exponential  rate  of  convergence  is  obtained  by  this  feedback 
process.  Finally  in  table  (6.1)  the  efficiency  index  of  the  error  predictor 
0  =  is  given  for  all  optimal  meshes.  0  ranges  from  0.61  to  1.42  which 

shows  that  the  relative  error  is  reasonably  estimated  by  the  error  prediction. 

We  performed  the  same  computations  on  geometrical  meshes  of  the  type 
shown  in  Fig.  3.1  with  elements  with  only  straight  edges.  Qualitative  and 
quantitative  behaviour  of  error  prediction  and  actual  error  was  very  similar 
to  the  meshes  with  circular  refinement  around  the  singularity,  so  that  these 
results  are  not  reported  in  detail. 

6.2.  Problem  2:  The  Wrench 

Figure  6.5  shows  the  domain  of  computation  and  the  loads  of  the  second 
problem.  The  domain  has  9  reentrant  corners  (including  the  two  tips  of  the 
crack  interior  to  the  domain).  Constant  traction  is  applied  along  the  edge  CD 
and  a  symmetry  boundary  condition  is  imposed  along  AB.  Isotropic  material 
with  E  =  3E8  psi  and  v  =  .3  was  assumed. 

Figure  6.6  shows  the  basic  mesh  of  this  problem  as  it  was  constructed 
interactively  by  a  modification  of  a  mesh  generator  of  MODULEF  [20].  The 
basic  mesh  is  as  coarse  as  possible,  only  modelling  the  geometry  of  the 
problem  (part  a)  in  section  4)) 


In  the  next  step  of  the  analysis  (part  b)  the  decision  about  critical  and 
noncritical  elements  is  made.  Again,  each  element  is  critical  as  each  has  at 
least  one  reentrant  corner  as  node. 

Step  c)  extracts  critical  data  about  the  solution.  As  was  shown  in 
section  5,  the  stress  intensity  factors  Kj  ^  and  exponents  a  ^  of  the 
singular  functions  provide  means  to  predict  the  error  in  the  energy  norm  for 
all  combinations  of  (geometric)  mesh  and  polynomial  degree  p.  For  a  suffi¬ 
ciently  accurate  extraction  of  stress  intensity  factors  it  is  necessary  to 
have  a  mesh  such  that  at  most  one  singular  point  is  in  each  element.  There¬ 
fore  out  of  the  basic  mesh  an  elementary  mesh  is  constructed  automatically, 
Fig.  6.7.  Table  6.2  gives  the  stress  intensity  factors  K,  ,  ,  1=1,2  for 

J  f  1 

d=2.  As  reference  solution,  "exact"  stress  intensity  factors  are  shown,  which 
were  obtained  on  a  strongly  refined  mesh  with  p=5  and  2985  degrees  of  free¬ 
dom.  It  can  be  seen  that  all  stress  intensity  factors  are  within  a  relative 
error  of  about  15?  compared  to  the  largest  stress  intensity  factor  (K2  -j  = 
-.15715*1  EM  at  point  1).  The  strain  energy  on  the  elementary  mesh  for  p=1 
is  U1  =  0.0323163736  and  for  p=2  U2  =  0.0357676705.  Using  6  =  .75  in 
(5.12),  the  estimated  exact  energy  is  U  0.03673856,  which  yields  an 
estimated  error  in  the  energy  norm  of  1M.M?  for  the  elementary  mesh  with 
p=3  (U3  =  0.0359766585).  This  value  is  used  to  calibrate  the  error  predictor 

(see  5.10). 

Step  d)  predicts  the  error  for  p=2,3,M  and  5  for  various  combinations  and 
optimizes  mesh  and  polynomial  degree  fo”  a  sequence  of  Nq  as  described  in 
section  5.2.  The  result  of  this  optimization  is  shown  in  Fig.  6.8  and  table 
6.3  as  it  is  presented  to  the  user  in  step  e).  Each  point  in  Fig.  6.8  corre¬ 
sponds  to  a  certain  layer  distribution  given  in  table  6.3.  Again,  the  next 
si'ep  in  practice  would  be  the  decision  of  the  user  on  a  particular  mesh-degree 


point 

factors 

basic  mesh 

refined  mesh 

1 

KI 

0.359377  E3 

0.377353  E3 

Kn 

-0.155661  E4 

-0.1571  54  E4 

2 

KI 

0.551  988  E3 

0.574465  E3 

KII 

-0.992989  E2 

-0.109476  E3 

3 

KI 

-0.536166  E3 

-0.534692  E3 

KII 

-0.205327  E2 

-0.983101  El 

4 

KI 

0.299322  E3 

0.539827  E3 

KII 

-0.943211  E2 

-0.168955  E3 

5 

KI 

0.561447  E3 

0.580995  E3 

KII 

0.127882  E3 

0.128434  E3 

6 

KI 

0.358611  E3 

0.374829  E3 

KII 

0.156916  E4 

0.155843  E4 

7 

KI 

0.325866  E3 

0.311223  E3 

KII 

0.115326  E3 

0.753*121  E2 

8 

KI 

0.420292  E3 

0.4373*15  E3 

KII 

-0.295953  E2 

0.224532  E2 

9 

KI 

0.204060  E3 

0.234956  E3 

*11 

-0.21  9939  E3 

-0.21  8943  E3 

Element  layers  of  point 


Polynomial 
I  degree  p 


Total  number 
of  degrees  of 
freedom 

1 

2 

3 

4 

5 

6 

7 

8 

9 

l 

488 

2 

2 

2 

1 

2 

2 

1 

1 

1 

692 

3 

2 

2 

2 

2 

3 

2 

2 

2 

861} 

3 

2 

3 

3 

3 

3 

2 

2 

3 

703 

2 

1 

2 

1 

1 

2 

1 

2 

1 

9111 

2 

2 

2 

2 

2 

2 

1 

1 

2 

1191 

3 

2 

2 

3 

2 

3 

2 

2 

2 

ino9 

3 

2 

3 

3 

3 

3 

2 

2 

3 

1673 

3 

3 

3 

4 

3 

3 

2 

3 

4 

21  87 

5 

1} 

}} 

4 

4 

5 

3 

3 

5 

1021 

2 

1 

1 

2 

1 

2 

1 

1 

1 

1229 

2 

1 

2 

2 

2 

2 

1 

1 

1 

1685 

2 

2 

2 

2 

2 

3 

2 

2 

2 

2101 

3 

2 

2 

3 

3 

3 

2 

2 

3 

2517 

3 

3 

1} 

4 

3 

3 

2 

2 

3 

3192 

4 

3 

4 

5 

3 

4 

3 

3 

5 

1127 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1955 

2 

1 

2 

2 

1 

2 

2 

1 

2 

2369 

3 

2 

2 

2 

2 

2 

1 

2 

2 

2985 

2 

2 

3 

3 

2 

3 

2 

2 

3 

Table  6.3.  Distribution  of  layers 
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combination.  In  contrast  to  that,  we  performed  here  finite  element  computa¬ 
tions  for  all  meshes  and  degrees  as  specified  in  table  6.3.  The  refined 
meshes  were  constructed  automatically  by  the  mesh  generator.  Figure  6.9  shows 


one  example  of  a  refined  mesh,  corresponding  to  the  combination  with  p=5  and 
2985  degrees  of  freedom  in  table  6.3.  Figure  6.10  shows  the  real  error  for 
all  computations.  We  estimated  an  exact  energy  from  extrapolation  on  an  ex¬ 
tremely  refined  mesh  and  high  polynomial  degree  to  UEX  =  0.0367809,  which 

is  very  close  to  the  estimated  energy  after  the  crude  computation  in  step  c). 
Comparing  Figs.  6.8  and  6.10,  it  can  be  seen  that  the  real  error  for  a  certain 
p-degree  levels  off  earlier  as  the  predicted  error.  This  seems  to  be  due  to 
the  fact  that  the  error  predictor  neglects  the  error  in  "smooth"  parts  of  the 
solution  which  are  still  significant  for  low  p-degrees. 

However,  8  of  10  optimal  points  (points  at  the  lower  left  envelope  in  Fig. 
6.8  which  are  marked  by  circles  in  Fig.  6.10)  turn  out  to  be  really  the  best 
meshes  for  the  specified  number  of  degrees  of  freedom,  the  remaining  two 
meshes  having  only  mildly  larger  error  than  the  actually  best  combinations. 

It  should  again  be  mentioned  that  the  optimal  meshes  constructed  in  the  feed¬ 
back  process,  converge  with  exponential  rate,  which  shows  that  the  feedback 
really  yields  optimal  mesh  design.  Finally  table  (6.4)  gives  the  efficiency 
index  0  *  —  for  all  the  meshes  which  is  in  the  range  of  0.43  to  1.35. 

By  this  feedback  process  it  is  possible  to  construct  a  mesh  which  yields 
an  accuracy  in  the  energy  norm  of  about  2%  for  3000  degrees  of  freedom. 
Extrapolation  from  the  elementary  mesh  shows  that  this  accuracy  could  be 
obtained  with  a  quasiuniform  mesh  using  linear  elements  (assumed  convergence 
rate  6  *  .25  )  with  about  1 E7  degrees  of  freedom.  An  adaptively  constructed 
mesh,  using  the  h-version  with  linear  elements  (convergence  rate  6  =  .5  ) 
should  yield  this  accuracy  with  about  34,000  degrees  of  freedom.  About  the 


Figure  6.9.  Refined  mesh  for  p=5,  2985  ZfDOF. 
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same  amount  of  dofs  would  be  necessary  for  a  pure  p-version  on  the  elementary 
mesh.  From  these  estimations  it  can  be  seen  that  about  one  order  of  magnitude 
in  the  number  of  degrees  of  freedom  is  gained  compared  to  an  adaptive  h- 
version  or  a  pure  p-version.  So,  roughly  two  orders-of-magnitude  in  storage 
and  computational  time  is  gained  if  high  accuracy  solutions  with  an  error  in 
the  range  of  1  to  2%  are  desired. 

The  additional  cost  for  the  feedback  is  negligible.  On  an  Apollo  420  the 
overall  time  for  the  computation  on  the  elementary  mesh  (p=1,2,3)  together 
with  graphical  mesh  generation,  the  sequence  of  optimizations  of  section  5.2 
for  p=2, 3, 4,5,  Ng  from  400  to  4000  in  steps  of  200  degrees  of  freedom  and  the 
construction  of  a  refined  mesh  with  2985  dofs  was  less  than  800  CPU  sec, 
whereas  the  final  computation  on  the  refined  mesh  with  p=5  took  more  than  4000 
CPU  seconds. 

Moreover,  it  should  be  mentioned  that  the  human  time  is  completely  inde¬ 
pendent  of  the  desired  accuracy,  a  situation  which  is  totally  different  from 
conventional  finite  element  analysis,  where  higher  accuracy  can  only  be 
achieved  by  time  consuming  construction  of  refined  meshes.  Only  the  basic 
mesh  has  to  be  designed  by  human  interference  and,  after  this,  only  decisions 
about  the  progress  of  the  analysis  have  to  be  made. 
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Polynomial 
degrees  p 

Total  number 
of  degrees  of 
freedom 

P(S) 

. 

me<*> 

e  1 

1 

2 

488 

17.21 

1  2.74 

1  .35 

2 

692 

9.95 

9.42 

1  .06 

2 

864 

8.07 

8.56 

0.94 

3 

703 

9.99 

12.54 

0.80 

3 

941 

5.58 

7.71 

0.72 

3 

1191 

4.32 

7.05 

0.61 

3 

1409 

3.34 

6.38 

0.52 

3 

1673 

2.92 

6.18 

0.47 

3 

21  87 

2.87 

6.17 

0.47 

4 

1021 

7.01 

8.83 

0.79 

4 

1229 

5.50 

7.46 

0.74 

4 

1685 

2.86 

4.39 

0.65 

4 

2101 

1.85 

2.99 

0.62 

4 

2517 

1  .26 

2.53 

0.50 

4 

3192 

1  .00 

2.30 

0.43 

5 

1127 

7.00 

9.33 

0.75 

5 

1  955 

3.67 

4.71 

0.78 

5 

2369 

2.24 

3.30 

0.68 

5 

2985 

1  .29 

1.94 

0.66 

Table  6.4.  Estimated  error,  exact  error  and  effectivity  index. 
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7.  CONCLUSIONS 


The  main  ideas  of  a  knowledge-based  system  advising  the  analyst  how  to 
design  the  mesh  and  degree  distribution  for  the  hp-version  of  the  finite 
element  method  have  been  presented.  The  system  includes  preliminary  analysis 
of  the  problems  and,  based  on  this  analysis,  it  advises  the  analyst  how  to  get 
the  prescribed  accuracy  for  the  lowest  cost.  Although  the  ideas  are  re¬ 
stricted  to  a  particular  class  of  problem  one  can  expect  that  ideas  of  the 
kind  we  presented  can  be  extended  to  more  general  cases  and  combined  with 
other  CAD  tools,  solid  modellers,  etc.,  to  a  larger  scale  expert  system. 
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The  Laboratory  for  Numerical  analysis  is  an  integral  part  of  the 
Institute  for  Physical  Science  and  Technology  of  the  University  of  Maryland, 
under  the  general  administration  of  the  Director,  Institute  for  Physical 
Science  and  Technology.  It  has  the  following  goals: 

o  To  conduct  research  in  the  mathematical  theory  and  computational 

implementation  of  numerical  analysis  and  related  topics,  with  emphasis 
on  the  numerical  treatment  of  linear  and  nonlinear  differential  equa¬ 
tions  and  problems  in  linear  and  nonlinear  algebra. 

o  To  help  bridge  gaps  between  computational  directions  in  engineering, 
physics,  etc.,  and  those  in  the  mathematical  community. 

o  To  provide  a  limited  consulting  service  in  all  areas  of  numerical 
mathematics  to  the  University  as  a  whole,  and  also  to  government 
agencies  and  industries  in  the  State  of  Maryland  and  the  Washington 
Metropolitan  area. 

o  To  assist  with  the  education  of  numerical  analysts,  especially  at  the 
postdoctoral  level,  in  conjunction  with  the  Interdisciplinary  Applied 
Mathematics  Program  and  the  programs  of  the  Mathematics  and  Computer 
Science  Departments.  This  includes  active  collaboration  with  govern¬ 
ment  agencies  such  as  the  National  Bureau  of  Standards. 

o  To  be  an  international  center  of  study  and  research  for  foreign 

students  in  numerical  mathematics  who  are  supported  by  foreign  govern¬ 
ments  or  exchange  agencies  (Fulbright,  etc.) 

Further  information  may  be  obtained  from  Professor  I.  Babuska,  Chairman, 
Laboratory  for  Numerical  Analysis,  Institute  for  Physical  Science  and 
Technology,  University  of  Maryland,  College  Park,  Maryland  207^2. 


